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ABSTRACT 

The  deterministic  Duane  Model  as  well  as  several  stochastic 
models  for  reliability  growth  are  studied.  Graphical  and 
quantitative  techniques  for  assessing  the  performance  of  these 
models  are  introduced  and  are  illustrated  on  a set  of  data.  A 
small-scale  simulation  is  carried  out  in  which  the  effect  of  two 
factors  (the  number  of  failure  modes  and  the  range  of  failure 
intensities)  are  investigated.  Some  practical  conclusions  are 
drawn  and  the  possibility  of  employing  doubly  stochastic  Poisson 
Processes  is  discussed. 


Key  Words:  Reliability  growth 

Piecevnse  constant  Poisson  process 
Nonhomogeneous  Poisson  process 
Duane  Model 


1 . I NTR05UCTI0N  AND  SUMMARY 


In  this  paper  we  carry  out  a comparative  study  of  various 
inodels  for  rel lability  growth  data.  We  assusss  the  reader  !s 
generally  familiar  with  the  problem.  For  a review  of  current 
work  in  th®  area  see  the  AKSAA  Rellsbillty  growth  Symposium  [19?43 

Of  these  snodels»  one  Is  deterministic  and  the  rest  stochastic 
The  former  Is  due  to  J.T.  Duane  [19643  and  is  widely  used 
In  practice.  It  seemed  of  some  interest  to  probe  the  relative 
strengths  and  weaknesses  of  the  two  approaches  and.  Indeed, 
Interesting  differences  were  found. 

Section  2 describes  the  models  and  discusses  their  motlvotlon 
In  Section  3 these  models  are  applied  to  a sot  of  reliability 
growth  data  and  various  goodness>of>f1t  procedures,  both  graphical 
and  quantitative,  are  introduced  and  applied.  Section  4 provides 
the  details  of  a Honte  Carlo  study  that  was  carried  out  to  shed 
more  light  on  how  well  the  different  models  performed  under 
various  conditions.  Section  5 analyzes  the  results  of  the  Konte 
Carlo,  while  some  tentative  conclusions  and  prospects  for  further 
study  are  given  in  Section  6. 

2.  THE  MODELS 

He  first  require  some  notation.  If  a device  has  been  on 
test  for  time  t and  has  experienced  H(t)  failures  up  to  that 
time,  then  Y(t)»t/N(t)  Is  called  the  cumulative  mean  time  between 
failure  (CRTBF)  at  time  t . If  the  failures  occurred  at  times 
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2. 

0«t^  < < tg  <...<  tyj  “ T,  then  we  define 

(f«l»2,....n}  as  the  n intsrfailure  times.  At  time  point 

t • EX|  can  be  interpreted  as  the  instantaneoL*s  mean  time 
1-1  ' 

between  failures  (IMTBF).  Contract  specifications  are  often 
phrased  in  terms  of  the  CHTBF  reaching  some  assigned  level. 
Nonetheless,  the  IMTBF  is  rather  important  because  It  gives  the 
rellebilltj^  engineer  a good  idea  of  current  progress  and  how 
long  until  the  CHTBF  goal  is  reached.  This  brings  us  to  the 
crulcal  point:  The  success  of  the  models  roust  be  judged  on  how 
well  they  allow  prediction  of  future  failures. 

The  deterministic  model  was  first  Introduced  by  J.T.  Duane 
Ci9643.  He  found  empirically  that  plots  of  CMTBF  versus  operat- 
ing time  were  approximately  linear  on  a log-log  scale.  That  Is, 

In  Y(t)  -v  a+b  In  t.  (1) 

Reliability  growth  occurs  when  Y{t)  Increases  with  t,  so  that  the 
slope  b is  positive.  To  calculate  the  IHTBF,  one  computes  the 
current  "intensity”  and  takes  reciprocals: 

I - (l-b)  e-®  f".  (2) 


“ dN(t)  'X*  d e"®  t^“^ 
'St 


1 


3 


or 

IMTBF  e*  . 


C3) 


Of  course  these  expressions  ought  not  be  given  a stochastic 
Interpretation.  Duane  suggested  using  least  squares  to  estimate 
a and  b In  (I).  Despite  the  obvious  theoretical  objections 
(after  all*  the  successive  observations  on  Y(t}  are  very  much 
dependent!)  the  procedure  has  worked  well  enough  In  practice  to 
achieve  wide  success. 

A stochastic  model  for  reliability  growth  was  proposed  by 
L.H.  Crow  [19753.  Using  (2)  as  a guide,  he  suggested  modelling 
the  process  as  a nonhomogeneous  Poisson  process  (NHPP)  with  an 
Intensity  function  of  the  form 
X(t)  • 

where  reliability  growth  Is  said  to  occur  If  B<1.  Naximuia  like- 
lihood estimates  of  y and  0 were  derived  as  well  as  formulas  for 
CMTBF  and  INTBF: 

EY{t)  • Y**  f®  (5) 


and 


1-8 

EX^  » (yD)'^  {1«2,3,...,n) 


EX 


1 


r(l/0+l). 


4. 

Note  that  since  the  Intensity  function  varies  with  time,  the 
successive  Interfailure  times  are  not  Independent. 

A second  stochastic  model • proposed  by  Braun  [19763  tool;: 
a different  tack.  One  can  think  of  reexpressing  the  NHS  of  (2) 
In  terms  of  N(t)  rather  than  t.  obtaining  something  of  the  form 

1-1/e 

X(t)  « re  CN(t)3  , (7) 


and  reliability  growth  occurs  If  e<l.  Thus  we  think  of  the 
failures  generated  by  a Poisson  Process  whose  Intensity  function 
changes  only  with  the  number  of  failures,  rather  than  with  time. 
As  a result  the  Interfailure  times  are  Independent  and 
exponentially  distributed  with  parameters 

1-1/6 

® X(t^)  “ y6  1 (1"1 ,2 » • • • tO ) . (8) 

Of  course  then 


EX 


and 


EY(t^)  «■ 


1 

I _±_ 


(9) 


(10) 


A variant  of  this  model  Is  obtained  by  changing  (8)  so  that  has 
a different  functional  dependence  on  1: 


I 


5. 

(11) 


X 


i 


c exp 


reliability  growth  occuring  If  b<0. 

This  has  an  advantage  of  allowing  the  Intensities  to  tend  to  a 
non>zero  constant»  a proposal  which  seems  quite  reasonable.  The 
parameters  of  these  models  can  be  easily  estimated  by  maxlmua} 
likelihood. 

The  estimation  theory  for  the  model (4)  Is  given  In  Crow  C13703, 
while  that  for  models  (7)  and  (11)  can  be  found  In  Appendix  A below. 

The  four  models  described  above  will  be  denoted  by  D«C»B1»BS 
respectively.  The  latter  three  despite  their  somewhat  ad  hoc 
parametrizatlons  seem  reasonable  competitors  of  model  0»  though 
they  In  no  way  attempt  to  model  the  physical  processes  generating 
the  failures.  Me  shall  see  In  later  sections  how  successfully  they 
accomplish  their  task.  Another  approach  to  the  problem  would  be  to 
use  techniques  analogous  to  Isotonic  regression.  This  would 
Involve  modelling  the  Intensity  function  as  a step  function  con- 
stant between  failures  and  required  to  be  non-increasing.  This  Is 
discussed  In  a report  by  Donelson  E1975].  It  may  well  be  that 
this  non-parametric  approach  could  then  be  followed  up  by  a para- 
metric one»  but  this  was  not  tried  here.  It  should  be  noted  that 
the  question  of  whether  to  use  N{t)  or  t as  the  Independent 
variable  In  stochastic  modelling  of  this  kind  Is  not  new.  See 
for  example  the  comments  of  T.  Lewis  In  the  discussion  of  a paper 
by  D.R.  Cox  [19553. 


6 


3-  ANALYZING  THE  DATA 

The  data  set  to  be  analyzed  here  consists  of  the  recorded 
failure  times,  during  development  testing,  of  a complex  electronic 
system.  The  data,  which  Is  presented  In  Appendix  B,  Is  composed 
of  52  failures  which  are  Indentifled  either  as  having  been 
generated  by  one  of  12  system  modes  (or  components)  or  as  a so- 

called  non-pattern  failure.  This  set  Is  denoted  by  SEl.  A 
subset  of  6E1,  called  6E2,  can  be  derived  by  only  recording  the 
first  failure  attributable  to  a particular  failure  mode  plus  all 
the  non-pattern  failures.  6E2  has  27  failures.  Table  I has  the 
parameter  estimates  for  all  the  models  for  both  data  sets.  One 
can  see  that  there  Is  agreement  among  the  models  that  6E2  dis- 
plays the  greater  reliability  growth. 

We  now  proceed  to  consider  some  descriptive  goodness-of-fit 
techniques  that  can  be  used  to  compare  the  suitability  of  the 
models,  with  the  exception  of  Model  D.  They  are  all  based  on  a 
probability  Integral  transform  of  the  Interfallure  times.  If 
has  cdf  F^,  then  Is  distributed . as  a uniform  (0,1) 

deviate.  In  our  case,  Involves  certain  parameters  which  must 
be  estimated  from  the  data.  As  a result,  the  collection 

|p^(X^):  i»l,2,...,n|  will  be  only  approximately  uniform  and 
Independent.  Nonetheless  we  can  try  plotting  F^(X^)  against  1 
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* 


to  see  whether  there  is  any  pattern  or  whether  the  are 

fairly  well  scattered  about  the  line  at  ,5,  (This  was  suggested 
by  J.H.  Tukey).  It  is  also  worthwhile  having  a QrQ.  plot  of  the 

A 

ordered  F^(X^)  against  i/n+1,  though  this  should  tend  to  look 
"super-regular",  A quantitative  feeling  for  the  informatsori 
contained  in  such  a graph  can  be  obtained  by  the  following 
strategy. 

A 

If  we  let  (i=l,2*...,n)  and  compare  the 

empirical  distibution  of  the  ^Z^|  with  that  of  the  uniform  dis- 

tributiont  then  it  Is  well  know  (cf  Durbin  C19733  for  example) 
that  the  usual  goodness-of-fi t statistics  such  as  the  Kolmogorov- 
Smirnov  do  not  have  the  same  distribution,  even  asymptotically, 

A 

that  they  do  when  the  true  is  used  in  place  of  F^.  How  it 
seems  plausible  that  if  we  were  to  randomly  partition  the  data  in 
to  groups  each  containing  say  10%  of  the  points,  and  carry  out  a 
goodness-of-fi t test  separately  for  each  group,  then  the  effect 
of  having  estimated  parameters  from  the  whole  data  set  should  be 
small,  when  considered  from  the  point  of  view  of  a single  group. 

In  fact,  we  can  use  the  methods  of  simulataneous  inference 
to  obtain  an  approximate  overall  level  of  significance  for  this 
procedure.  A theoretical  justification  of  this  approach  Is 
currently  under  preparation  by  the  first  author. 
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Crow  C19753  has  adapted  the  work  of  Darling  C19SS]  and  shown 
how  to  make  use  of  the  Cramler-von  Kises  statistic  for  testing  the 

adequacy  of  Model  C with  parameters  estimated.  He  has  shown  that 
If  the  model  Is  correct*  then  the  distribution  of  the  statistic  Is 
Independent  of  the  actual  values  assumed  by  the  parameters  and  has 
by  Monte  Carlo  computed  the  critical  values  of  the  statistic.  Note 
that  this  approach  works  only  for  Model  C*  whereas  the  method 
suggested  previously  Is  quite  general  In  Its  applicability. 

A 

Figure  1 contains  a plot  of  F(X|}  versus  1 for  Model  BI  with 
GEl  data  and  Figure  2 contains  a plot  of  F(X^)  versus  1 for  • 
Model  C with  GE2  data.  These  plots  as  well  as  others  not  shown 
here  are  reasonably  patternless*  having  a hint  of  a positive  slope 
which  would  suggest  that  the  (stochastic)  models  tend  to  over- 
estimate the  earlier  failure  times  and  underestimate  the  later  . 
ones.  The  QrO.  plots  are  all  well  behaved.  The  modified  Kolmpgorov- 
Smlrnov  statistic  was  computed  for  each  case  and  was  never  found 
significant  (o».05).  In  addition,  the  modified  Cramler-von  Mises 
statistics  was  calculated  for  Model  C and  was  not  significant 
(o".10).  Finally,  Figure  3 plots  In  y(t)  vs  In  t for  the  GEl  data. 
Note  that  It  Is  a bit  u-shaped.  Indicating  that  CMTBF  does  not  In- 
crease until  the  process  Is  well  under  way. 

What  have  we  learned?  Only  that  both  data  sets  exhibit  re- 
liability growth*  6E2  more  so  than  GEl,  and  that  the  models  seem 
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to  do  an  equally  good  job  of  fitting  the  data.  In  fact,  our 
graphical  and  quantitative  techniques  so  far  give  us  little  reason 
to  choose  between  the  models.  It  would  be  worthwhile  to  see 
whether  the  forms  of  the  different  models  have  consequences  which 
can  be  used  in  an  exploratory  (graphical)  way. 

Suppose  we  let  denote  the  expected  value  of  , the 

interfailure  time.  Then  under  Model  Bl, 


1-1/B 


or 

Under  Model  C,  on  the  other  hand. 


(X2) 


1/pT-l) 


So  approximately. 


(13) 


Of  course,  (13)  can  be  obtained  approximately  from  (12)  by  expan 
Sion  of  the  logarithm.  Nonetheless  it  may  be  that  one  approach 


% 


.p 
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is  better  In  dealing  with  real  data.  Consider  now  replacing 
expected  values  in  (12)  and  (13)  by  the  actual  observations.  If 


Model  B1  holds,  then  plotting  log  against  1/1  should 

yield  a linear  regression  through  the  origin  with  slope  -(l-l/S). 
In  fact  the  estimated  slope  could  be  used  to  obtain  an  estimate  of 
3.  If  Model  C holds,  then  plotting  against  1/1  should 


yield  a linear  regression.  This  plan  was  carried  ou.t  with  in- 
teresting, but  by  no  means  conclusive  results. 

Using  the  GEl  data  as  a base,  the  Interfailure  times  /xA  were 


smoothed  and  from  the  smoothed  sequence  |X^,|the  sequences 

- Xi.i/xi  Md  |v,  - log  (x,.,/x,)j  constructed.  These 


were  in  turn  smoothed  (to  reduce  negative  correlations  between 

/-I  f I 

successive  elements)  yielding  jU.||  and  ]V.|j. 

I 

Regression  of  on  1/i  Indicated  a negative  slope  while 


regression  of  iV^fon  1/i  gave  a slope  estimate  which  when  solved 


for  3 gave  a value  of  0.72,  quite  close  to  the  maximum  likelihood 
estimate.  However  in  neither  case  was  the  regression  significant, 
so  the  evidence  in  favor  of  Model  B2  is  by  no  means  overwhelming. 
Nonetheless,  it  suggests  that  this  technique,  properly  enhanced 
by  smoothing,  may  be  of  some  value. 
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We  now  Introduce  three  quantities  which  attempt  to  measure 
how  well  the  models  do  at  predicting  both  the  IMTBF  and  the 
CMTBF.  These  measures  will  also  allow  us  to  compare  the  determi- 
nistic model  with  the  stochastic  ones  on  a somewhat  equal  footing. 
Now  stolctly  speaking,  we  should  reestimate  the  parameters  of  the 
model  after  each  failure  and  use  these  new  estimates  to  form  our 
prediction  of  the  next  failure  time.  This  would  be  rather  time 
consuming,  so  Instead  parameters  estimated  from  the  entire  data 
set  were  used  to  make  the  predictions.  However  for  the  6E1  data 
the  scheme  of  successive  reestimation  was  carried  out  and  the 
parameter  estimates  remained  fairly  stable,  though  less  so  for 
Model  B2. 

The  statistics  described  below  are  closely  related  to  the 
criteria  used  In  the  study  described  In  Schafer,  etal  C1975J. 

t / (n-2) 

1«1  ’ ’ 


” (X.-X)^  / (n-1) 

1«1  ^ 


_ n 

X « E X./n 
1“1  ’ 


The  first  is 


R1  - 


where 
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and 


X| 

'o 


(x) 


The  motivation  for  this  Is  obvious,  and  we  expect  R1  to  be  small 
(certainly  less  than  one)  If  the  model  Is  of  any  use.  A variant 
of  Rl,  denoted  by  R2,  Is  meant  to  reduce  somewhat  the  effect  of 
the  great  Irregularity  In  successive  Interfallure  times  on  Rl* 

We  begin  at  some  appropriate  time  once  the  process  Is  underway 
(say  after  10  failures)  and  collect  the  failure  times  In  successive 
groups  of  those  after  that  point.  Let 

m = number  of  groups  of  three 

'V. 

Xj  • mean  of  the  observed  for  the  j group 

EXj  “ EX^  for  the  middle  time  In  the  group 


•V 

X = 


and  define 
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in  *'*9 

2 (X.-  EX.)^  /(m-a) 

R2  j»i  3 J 

m 'u  - 

Z (X.-X)^  /(m-l) 


The  third  criterion,  used  In  Schafer,  etal  [19753  deals  with 
prediction  of  CMTBF.  Let  Y(t^)  denote  the  estimate  of  Y(t|),  the 

CMTBF  at  time  t^ . Then 


£ lY(t,)-Y(t.y  /(n-2) 
R3  - f-1  ’ ’ 


£ (Y(t,)-T(t,))^  /(n-1) 


n 

where  ^(t^)  * I Y(t.)/n. 

1 1-1  ’ 

Table  II  contains  the  results  of  computing  these  statistics 
for  each  of  the  models  and  data  sets.  It  seems  clear  that  the 
stochastic  models  do  better  with  the  IMTBF  while  the  deterministic 
one  does  better  with  the  CMTBF,  but  there  Is  little  to  choose 
between  the  stochastic  models.  A few  general  points  may  be  worth 
mentioning.  Since  CMTBF  Is  more  stable,  all  the  models  do  well 
with  respect  to  the  R3  criterion,  and  better  with  the  GE2  than 
GEl  data.  This  suggests  that  looking  at  only  first  Instances  of 
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failure  inodes  can  be  one  way  of  obtaining  better  predictions. 

Of  course  there  is  a practical  difficulty  that  failures  can  not 
always  be  correctly  classified  or  one  not  even  recognized  as 
pattern  failures  until  much  later.  Secondly*  a comparison  of  the 
R2  values  with  R1  values  shows  that  once  we  have  smoothed  some  of 
the  roughness  away,  none  of  the  models  does  better  than  the  sample 
mean.  This  is  most  pronounced  with  the  6E2  data  and  the  1n> 
escapable  conculslon  Is  that  much  more  Insight  Is  required  before 
we  can  do  better.  To  this  end,  a simulation  was  performed  which 
aimed  to  discover  what  aspects  of  the  process  had  the  greatest 
effect  on  the  performance  of  the  models.  This  Is  described  in 
the  following  section. 

4.  THE  SIMULATION 

The  simulation  was  designed  to  produce  data  that  could  be 
reasonably  considered  to  approximate  reliability  growth  data. 

An  observed  process  of  this  kind  results  from  the  superposition 
of  the  failure  patterns  of  different  modes  (usually  due  to  design 
faults)  together  with  a process  of  non  pattern  failures  (wear  outs, 
breakages,  etc.)  that  Inevitably  occurs  during  the  operating  life 
of  any  system.  Since  the  6E1  data  was  used  as  a guideline,  let 

us  review  Its  structure.  There  are  thirteen  failure  modes  which 
when  treated  separately  as  homogeneous  polsson  processes  yield 
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estimated  Intensities  of  between  .5  and  2.5  failures  per  1000  hrs. 
(f/th).  The  non  pattern  failures  have  an  estimated  intensity  of 
4.4  f/th.  The  number  of  observed  failures  per  mode  varies  between 
two  and  seven*  with  most  being  either  two  or  three. 

Our  simulated  process  will  then  be  constructed  as  the  super- 
position of  a number  of  independent  Poisson  processes*  each  to  be 
truncated  after  a random  number  of  failures.  Even  so  simple  a 
process  defies  analytic  Investigation.  It  resembles  the  so-called 
"branching  Poisson  process"  Introduced  by  Lewis  [19643  as  a model 
for  computer  system  failures.  His  process  was  stationary  but  even 
so,  the  main  results  were  all  asymptotic.  For  our  case  we  are 
precisely  Interested  In  the  transient  phase  and  are  not  concerned 
with  the  structure  of  the  non-pattern  failures  so  asymptotics 
would  by  meaningless  here.  Nonetheless,  the  work  on  superposition 
of  Independent  sparse  processes  (see  the  review  by  Cinlar  [19723  ) 
makes  It  plausible  to  assume  that  the  superposition  will  resemble 
a non-homogeneous  Poisson  process. 

The  advantage  of  this  approach  Is  that  the  generating  scheme 
Is  based  on  an  analysis  of  the  physical  situation  and  not  on  a 
preconceived  model.  Hence  if  one  of  the  models  does  very  well, 
we  can  have  a measure  of  confidence  In  Its  success.  Of  course 
with  a Monte  Carlo  we  can  Investigate  certain  problems  which 
would  be  difficult  if  not  impossible  to  do  with  only  real  data. 
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The  assumption  of  Independence  In  the  model  can  not  be 
strictly  justified  In  practice  but  Is  probably  not  a bad  approxi- 
mation to  reality.  We  are  primarily  Interested  In  the  effect  of 
two  factors:  the  number  of  failure  modes  and  the  range  of  failure 
Intensities.  In  our  simple  two-way  design  each  factor  has  three 
levels.  For  the  number  of  failure  modes:  low  (5)>  medium  (XO), 
high  (15).  The  Intensity  ranges  are  low,  medium,  and  high,  chosen 
In  the  following  way: 

Let  m denote  the  number  of  failure  modes  (5,10,  or  IS). 

Let  Z|<Z2<. . denote  the  expected  values  of  the  order  statistics 
from  a uniform  (0,1)  random  sample  of  size  m.  The  intensities 

^^1^1»1  chosen  so  that: 


(a) 

low: 

>1  * 

2Z^+.5 

(1°1,2,. ..,m} 

(b) 

medium 

"l  “ 

2yTJ+.5 

(1=1,2, .,.,m) 

(c) 

hi  gh 

2 yz^+.5 

(1=1,2, . . . ,m). 

(There  Is  clearly  overlap  among  the  ranges).  We  then  generate 
ro  Independent  observations  ni***»*^n,  from  a Poisson  distribution 
with  parameter  three  using  the  generator  described  by  Ahrens 
and  Dieter  C197AT.  The  n.  are  randomly  assigned  to  the  . 
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Thus  we  generate  m Independent  Poisson  processes,  where 
th 

the  1 process  has  Intensity  X^.  and  Is  truncated  after  the 

appearance  of  the  n^^^  failure.  These  are  then  superposed  and 
to  them  Is  added  another  Independent  Poisson  process  with  an 
Intensity  of  3 f/th,  representing  the  non  pattern  failures.  The 
result  Is  a simulated  reliability  growth  process  which  Is  trun- 
cated for  fitting  and  analysis  at  the  first  failure  after  4500 
hours.  For  each  combination  of  the  two  factors  we  have  20  re- 
plications, so  that  there  are  3x3x20*180  failure  processes  genera- 
ted and  analyzed. 

It  should  be  emphasized  that  the  twenty  replications  In  a 
particular  cell  have  the  same  number  of  failure  modes  n . the 
1^^  mode  having  the  same  intensity  X^.  However  the  number  of 

A K 

failures  n^  observed  In  the  1 mode  Is  generated  anew  for  each 
replication.  Of  course  the  waiting  times  themselves  are  generated 
by  independent  exponential  variates.  The  simulation  was  pro- 
grammed In  FORTRAN  on  a POP  11/40  under  the  UNIX  operating  system. 
A (hopefully I)  trustworthy  congruentlal  generator  was  used  to 
produce  the  uniform  deviates  employed. 
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5‘  ANALYSIS  OF  RESULTS 

Three  models  were  fitted  to  each  replication:  0,C»  and  B2 
and  the  three  criteria  Rl,  R2,  R3  were  calculated  for  each.  (For 
reasons  of  economy.  Model  B1  was  Jettisoned  at  this  stage.)  A 
discussion  and  Interpretation  of  the  resulting  analyses  of 
variance  forms  the  body  of  this  section.  Before  doing  this.  It 
seems  useful  to  mention  some  of  the  general  Impressions  gained  by 
use  of  the  graphical  techniques  Introduced  In  Section  3,  though  It 
would  be  obviously  unwieldy  to  go  Into  too  great  detail. 

First,  all  the  models  did  better  the  larger  the  number  of 
failure  modes  and  the  higher  the  range  of  Intensities  > clearly 
the  effect  of  having  more  data  to  work  with.  The  stochastic  models 
had  a tendency  to  underestimate  the  amount  of  reliability  growth, 
I.e.  overestimating  the  earlier  failure  times  and  underestimating 
the  later  ones. 

Second,  the  plots  of  In  Y(t)  against  In  t often  had  a U-shaped 
appearance.  Initially  decreasing  but  then  Increasing  once  the 
process  was  well  under  way.  This  suggests  an  Interesting  modi- 
fication of  Duane's  original  proposal;  namely,  estimating  the 
parameters  using  weighted  least  squares,  with  the  weights  assign- 
ed to  In  Y(t)  proportional  (say)  to  /N(t).  This  would  certainly 
give  proper  due  to  the  Increasing  number  of  observations  contri- 
buting to  successive  values  of  Y(t}.  This  proposal  was  not, 
however,  tried  In  the  present  study. 
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These  graphical  techniques  and  even  the  goodness-of>f1t 
statistics  discussed  In  Section  3 were  not  sufficiently  delicate 
to  discriminate  between  the  models  and  so  we  turn  to  the  more 
quantitative  methods.  For  each  criterion  (R1,R2,R3}  a three 
factor  ANOVA  was  carried  out,  employing  the  criterion  as  the 
dependent  variable  with  the  factors  being  models,  number  of  modes, 
range  of  Intensities. 

In  addition  to  the  usual  ANOVA  by  least  squares  fitting, 
a three-way  median  polish  was  carried  out  (cf..  Tukey  ri9763). 

The  results  were  qualitatively  the  same,  though  the  magnitudes  of 
the  fitted  effects  differed  somewhat;  there  was  little  Indication 
that  a transformation  of  the  data  was  called  for.  We  therefore 
content  ourselves  with  presenting  only  the  normal-theory  analysis, 
and  Table  III  contains  a compact  summary  of  the  results.  After 
considering  these,  we  will  discuss  some  of  the  more  Interesting 
two  way  tables.  The  rather  large  number  of  degrees  of  freedom 
(df)  for  Residual  makes  the  analyses  fairly  clearcut. 

For  Rl,  the  main  effects  are  all  highly  significant,  and  the 
roodes-frequency  Interaction  Is  very  significant.  Model  D does 
quite  a bit  worse  than  C or  B2.  The  differences  among  the  levels 
within  the  other  two  factors  are  In  the  expected  direction  but  one 
not  as  large  In  magnitude  as  for  the  first.  For  R2,  only  models 
and  modes  are  significant.  Note  that  the  residual  mean  square  Is 
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larger  than  that  for  R1  by  a factor  of  300.  Again  Model  D does 
substantially  worse.  For  R3,  only  models  Is  significant,  but  here 
Model  0 does  substantially  better  than  the  others  and  Model  C 
substantially  worse. 

As  we  are  led  to  expect  from  Section  3,  Model  0 does  best 
with  CMTBF,  while  Models  C and  B2  do  best  with  IHTBF.  It  Is 
Interesting  that  as  the  "smoothness”  of  the  criterion  considered 
Increases  {Rl->R2-*-R3) , the  Importance  of  the  factors,  modes  and 
frequencies,  diminishes.  Let  us  Investigate  this  further,  by 
considering  some  two-way  analyses. 

Tables  IV,  V,  and  VI  contain  the  ANOVAs  and  fitted  effects 
for  Rl,  R2  and  R3  respectively.  There  we  note  that  Model  D Is 
sensitive  to  the  number  of  modes  for  Rl  and  R2  but  not  R3.  Model 
C is  sensitive  to  the  number  of  modes  for  all  criteria,  but  to  the 
frequency  range  only  for  Rl.  Model  B Is  sensitive  to  the  number 
of  modes  for  Rl  and  R2  and  to  the  frequency  range  for  Rl, 

However  when  we  consider  the  data  by  criterion,  we  see  that 
Model  C tends  to  be  quite  a bit  more  sensitive  than  the  other 
models  to  the  number  of  modes,  while  Models  B and  0 arc  roughly 
equivalent  in  this  regard.  It  is  to  facilitate  this  comparison 
that  the  tables  are  arranged  as  they  are.  Finally  nowhere  are  the 
two  factor  interactions  of  any  Importance  except  possibly  with 
Model  C for  Rl. 
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6.  CONCLUSIONS 

There  are  several  conclusions  that  can  be  drawn  from  the 
study.  First  the  graphical  techniques  are  quite  useful  In  giving 
a general  impression  of  the  fit,  but  are  not  delicate  enough  to 
distinguish  between  models.  Second,  smoothing  the  data  may 
certainly  enhance  the  stability  of  the  estimation  procedures  and 
give  some  hope  of  developing  ad  hoc  approaches  to  model  develop* 
inent.  The  Duane  model  is  quite  workable  as  far  as  CMTBF  Is 
concerned  but  its  performance  would  probably  be  Improved  by  use  of 
weighted  least  squares.  It  is  not  recommended  when  IMTBF  Is  of 
Interest,  there  the  stochastic  models  are  superior.  In  addition, 
data  of  the  GE2  type  seem  to  allow  regression  parameters  of  the 
model  to  be  estimated  more  accurately  than  data  of  the  6E1  type. 

Turning  to  the  simulation,  it  seems  clear  that  although  the 
various  techniques  work  better  with  increasing  amounts  of  data,  a 
really  successful  method  would  be  fairly  insensitive  to  the  factors 
considered  such  as  number  of  modes  and  range  of  intensities.  CThis 
Is  particularly  true  in  the  case  of  R3  which  deals  with  a stable 
criterion,  namely  CMTBF.)  In  fact  it  is  a bit  surprising  that 
these  factors  do  not  have  a greater  influence.  Of  course  a larger 
range  of  levels  ought  to  be  used  in  any  further  simulation.  Over- 
all Models  C and  B2  are  roughly  equivalent  in  performance,  though 
Model  B2  seems  less  sensitive  to  changes  In  the  levels  of  the 
factors. 
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We  are  now  In  the  situation  of  having  the  problem  somewhat 

In  hand,  but  by  no  means  mastered.  One  approach  would  be  to 

consider  stochastic  models  with  different  parametlzatlons.  The 
astute  reader  will  have  noticed  that  the  models  discussed  here 
bear  more  than  a passing  resemblance  to  the  linear  logistic  models 
of  Cox  [1970],  where  It  Is  the  log  odds  ratio  that  Is  posited  to 
have  a linear  regression  on  one  or  more  Independent  variables.  In 
the  present  context  It  seems  easier  to  work  with  the  discrete 
N(t)  rather  than  the  continuous  t Itself.  (See  Brown  [19723  for 
a discussion  of  hypothesis  testing  for  nonhomogeneous  Poisson 
processes. ) 

Thus  we  can  envision  trying  models  of  the  form 

log  X|  ® a+b  log  1 + c(log  1)^  (14) 

trying  to  get  a "best  fit".  It  does  not  seem  likely  that  this 

would  result  In  substantial  gains.  For  If  one  compares  the 
sequence  of  predicted  Interfailure  times  with  the  actual  observed 
sequence.  It  Is  clear  that  the  former,  while  capturing  the  drift 
of  the  latter,  does  not  resemble  It  In  the  least  - particularly 
because  of  the  great  variability  In  the  observed  times.  This 
situation  can  not  much  change  by  the  addition  of  terms  Into  the 
regression . 
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these  considerations  lead  us  to  suggest  that  perhaps  a doubly 
stochastic  process  model  might  be  more  suitable  (see  Srandell 
[19723.)  That  Is. 


log  ■ a+b  log  1+e^ 


(15) 


where  n[o,  a^J  and  Independent  for  each  1,  Such  a model 

would  certainly  Introduce  more  variability,  but  the  question 
arises  of  whether  we  understand  something  better  for  having  named 
It.  The  answer  In  this  case  Is  yes,  because  such  a model  has 
ramifications  which  can  be  sought  In  data.  If  Is  exponentially 
distributed  with  parameter  then  under  (14) 


(X4)  » 1 

^1 


while  under  (15) 


Var  (X^)  « E 


Var  11 

I-T7 


'w  1 + Var  I 

T?  1 X 

^1 


1 


Z4 


Thus,  In  this  latter  case,  there  Is  another  cotcponent  to  the 
variance. 

Unfortunately,  In  any  single  process  we  have  only  one 
observation  on  X^.  But  It  Is  quite  common  to  test  three  or  four 
devices  simulataneously,  thus  giving  replications  which  allow  us 
to  compute  estimated  variances.  Note  that  this  Is  only  possible 
when  we  use  N{t}  as  the  Independent  variable  (as  In  Models  B1  and 
B2}  rather  than  t ( as  In  Model  C). 

. For  definiteness,  consider  making  m replicate  observations 
pi^^M * particular  reliability  growth  process.  We  then 

have  n (say)  sets  of  lid  observations  fu  '■I?.! 

where  'x*  EXP  (X^)  . . . ,m. 

2 2 

Then  compute  Jj  • 1 Ex.,  and  s.  » 1 E . If 

’ "a  J ’ 


we  fit  a model  like  (7)  to  the  data  (using  the  full  likelihood). 


then  the  sequence 


W ■ 

I ^il 


w 


should  probably  conform  very  closely  to 


But  the  sequence 


would  probably  not  resemble 


very  much.  If  such  were  the  case  It  would  be  strong 


evidence  In  favor  of  (15),  particularly  If  the  regression  of 
2 

s<  - 1 on  1 had  zero  slope.  If  this  regression  were  not 

1 A ai 
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constant.  It  would  of  course  Indicate  modifying  (15)  to  allow 

2 

the  E^'s  to  have  different  variances  o^  , and  then  one  might  try 

2 

to  study  the  structure  of  the  cr^  . But  in  this  situation  there 
would  be  also  ample  Incentive  to  search  for  another  regression 
which  would  do  as  well  in  fitting  the  while  keeping  the 
2 

0^  constant. 

It  is  certainly  premature  to  proceed  further  along  these 
lines  until  supporting  evidence  for  models  like  (15)  can  be  found. 
But  if  they  were  to  prove  plausible,  they  would  provide  reliabi- 
lity engineers  with  better  variability  estimates  with  which  to 
plan  development  programs  and  the  like. 
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APPENDIX  A 

For  Model  B1  we  have 


X(t)  » . 

The  likelihood  function  of  the  observed  Interfallure  times  Is  given 
by 


l.{x,  X ) » IT  X.  e *1  ’ 


1»1 


where 


X^  » yBI 


1-1/B 


taking  H(o)  ■ 1 for  convenience.  Differentiating  the  log  likeli- 
hood with  respect  to  y and  B and  setting  the  resulting  expressions 
equal  to  zero  yields 


Y 


n 


r BX.1 
1 ’ 


1-1/S 


(Al) 


and 


1-1/S 

Z 1 x^  log  1 Z log  1 

1 " 1 

^ n 

1-1/S 

Zx^l 
1 ’ 


(AZ) 


Equation  (AZ)  can  be  easily  solved  numerically  for  6,  the 

A 

solution  being  used  then  In  (Al)  to  obtain  y.  Note  that  various 
quantities  of  Interest  are  now  easily  computed.  For  example 
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E(X^^l)  Is  estimated  by 


1 


8 


» 


and  from  this  one  can  calculate  the  expected  number  of  failures 
that  will  occur  in  a given  amount  of  time.  In  planning  a develop- 
ment program,  one  may  want  an  estimate  of  the  number  of  failures 
that  will  be  needed  to  bring  the  reliability  to  a desired  level. 

If  the  desired  intensity  is  and  the  trial  values  of  the 
parameters  (perhaps  current  estimates)  are  *»>d  3^^  > then  the 
required  number  is  approximately 


Finally,  one  would  like  to  have  approximate  standard  errors 
for  the  parameter  estimates.  If  we  let 


9 = 


and  1 (6^)  « log  L(e^;  Xj 
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then 


3^  1(0) 


30! 


3^  Hi) 


30: 


3^  He)  * 


30^302 


-n 


~1  E log  1 + 

e^Ex^dog  1)^  1 


-E  X. 
1 ^ 


[1  + 1 log  ll 

*2  J 


(A3) 


Then  as  Is  well  known  from  maximum  likelihood  theory*  (even  though 
the  observations  are  not  Identically  distributed) 

\ 


hi  <a) 


^3^  Hi) 


30^  3ej 


(A4) 


/ 


which  Involves  replacing  x^  In  equations  (A3)  by  X,j 

1 -1 


-1 


®1  ®2 


1-1/00 
1 ^ 


. . The  resulting  expression  depends  on  the 


.A  A 


data  only  through  e^,  @2  ond  n,  and  Is  particularly  simple  to 

A 

compute.  The  asymptotic  variance-covariance  matrix  of  0 Is  the 
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inverse  matrix  to  (A4),  denoted  by 
a consistent  estimate,  (e))  i 

For  Model  82  we  have 

X,  - c exp  |b  (l- 


(e)) 


and  which  has  for 


Differentiating  the  log  likelihood  with  respect  to  b and  c 
and  setting  the  resulting  expressions  equal  to  zero  yields 


c « n/E  x^  exp 


{- 


and 


Ex^jexp  b 


(-  i?)' 


(-  ¥)‘ 


Ex^  exp  )b 


-J  L.  J 


Equation  (AS)  can  be  easily  solved  numerically  for  b, 
being  used  in  (A5)  to  obtain  c . 


Now  let  9 « / 1 « 

0, 


and 


(A5) 


(A6) 


1(e)  = log  L(^;  Then 


the  solution 


The  argument  after  equations  {A3)  applies  here  as  well,  and  a 
consistent  estimate  of  the  variance-covariance  matrix  of  8.  Is 
obtained.  In  this  case  the  variance-covariance  matrix  depends 

A 

only  on  6^,  and  n. 
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A.  6E1  Data.  Interfallure  Times  (to  be  read  across  rows) 


25 

15 

210 

25 

20 

5 

15 

30 

5 

25 

75 

10 

110 

10 

80 

120 

60 

110 

10 

60 

30 

25 

175 

175 

25 

200 

175 

25 

10 

65 

25 

250 

5 

95 

50 

25 

45 

5 

15 

60 

70 

10 

170 

20 

80 

30 

195 

125 

100 

150 

60  190 


GE2  Data,  Interfallure  Times 

(to  be 

read  across  rows) 

25(2) 

15(NP) 

210(NP) 

25(NP) 

20(3) 

5(2) 

15{NP) 

30(7) 

5(NP) 

25(2) 

85(NP) 

IIO(NP) 

lO(NP) 

260(2) 

180(2) 

55(NP) 

175(2) 

175(NP) 

500(1) 

25(3) 

255(NP) 

215(3) 

5(NP) 

15(5) 

130(NP) 

200(NP) 

740(NP) 

' 

NOTE: 

Quantities 

In  brackets  give 

the  number  of  failures 

that  occurred  In  that  mode  whose  first  occurrence 

Is  recorded  In  the  corresponding  entry.  "NP"  stands 
for  non-pattern  failure. 


TABLE  1:  PARA?1ETER  ESTIMATES 

Data 


Model  ^N.  GEl  GE2 


A 

A 

b - .195 

b - .347 

D 

a * 2.65 

a • 1,89 

A 

A 

B = .837 

B » .620 

C 

A 

A 

Y « .0537 

Y * -171 

A 

A 

B * .799  (.1002) 

B * .592  (.00817) 

Y « .0385  (.0236) 

Y * .0799  (.00577) 

b * .635  (.457) 

b - 1.49  (.0197) 

B2 

c « .0121  (.00253) 

— — 

c « .0004  (.000117) 

NOTE:  Quantities  in  brackets  are  estimated 

standard  errors  of  estimates. 
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TABLE  3:  3-FACTOR  ANOVA'S  AND  ESTIMATED  MAIN  EFFECTS  (FROM 


NOTES:  (1)  Significance  levels  of  F-statistics ; ♦-^p<.05; 

(2)  Bracketed  quantities  in  "Residual"  row  are  mean 
squares.  Thus  all  mean  squares  and  sums  of  squares 
can  be  computed  from  information  in  Table 

(3)  e.g.  Estimated  main  effect  of  B2  with  Rl  as  dependant 
variable  is  .92+(- .04)=.88 


TABLE  5:  2-FACTOR  ANOYA'S  WITH  R2  AS  DEPENDE?4T  VARIABLE;  ALSO 
ESTIMATED  MAIN  EFFECTS 


Grand  Meanj  I 
'5 

Modes  ■ 10  2 

15 


^Low 

I ? * »3» 


. 

High 


MODEL 


as  X 

Freq  4 , 18 

MS 

Residual  171  (7.02) 

Total  180 


Est'd 

Effect 


1.72 

.70 

3.17*  -.37 

-.33 
-.27 

.46  .10 


MODEL 


F-  Est' 

Stat  Effe< 


1.10 

.19 

7.04**  -.09 

-.10 
-.03 


H3 

(.234) 


MODEL 

B2 

F- 

Stat 

Est'd 

Effect 

1.28 

.40 

3.14 

-.15 

-.?5 

1 

.27 

1.44 

-.15 

-.12 

1.43 

MS 

(2.31) 

TABLE  6:  2 FACTOR  ANOVA'S  WITH  R3  AS  DEPEHDERT  VARIABLE 


Source 

df 

0 

1 

C 

B2 

Grand  Mean 

1 

- 

1 

- 

Modes 

2 

.24 

3.45* 

in 

o 

• 

Freq 

2 

1.66 

1.96 

.13 

Modes  X 

Freq 

4 

1.02 

1 

.54 

l.U 

MS 

4 

(K  j 

m 

1 

■371  * 
1 " - 1 

(.120)  j 

i / 

{.r/3} 

